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Abstract 

Two methods for computing bound states of the Helmholtz equation in a finite 
domain are presented. The methods are formulated in terms of the Dirichlet-to- 
Neumann (DtN) and Neumann-to-Dirichlet (NtD) surface integral operators. They 
are adapted from the DtN and NtD methods for bound states of the Schrodinger 
equation in R 3 . A variational principle that enables the usage of the operators is 
constructed. The variational principle allows the use of discontinuous (in values or 
derivatives) trial functions. A numerical example presenting the usefulness of the DtN 
and NtD methods is given. 

1 Introduction 

In this work we are interested in finding solutions to the problem consisting of the 
Helmholtz equation defined in a two- (or more) dimensional finite volume T: 

Atf(r) + £; 2 ^(r) = 0, r G T (1) 

and the homogeneous Dirichlet condition on the boundary dT: 

*(r) = 0, re dT. (2) 

The set ([T|- (|2|) is an eigenproblem in which the values {— A; 2 } are the eigenvalues and 
{^} are the corresponding eigenfunctions. The eigenproblem (pQ)-([2]) appears in many 
different areas of physics. It describes, for example, the behaviour of a particle confined 
in an infinitely deep potential (in this case k 2 is proportional to the energy of the particle 
while |^| 2 is the probability density) or vibrations of a homogeneous membrane (k is 
proportional to the vibration frequency, \E' is the amplitude), it is useful in studying the 
propagation of electromagnetic waves in waveguides, etc. So, although the problem of 
finding eigenvalues and eigenfunctions of the Laplace operator has been known for many 
decades, it remains very important in many fields. 

The standard analytical approach to problems like (HJ-© is the method of separation 
of variables. The first step in the method is to choose an appropriate coordinate system. 
The choice depends on the shape of dT. In practice, only in some cases it is possible to 
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find a system fitted to the geometry of a problem and to obtain the exact solutions using 
the separation of variables technique. In general, the shape of the boundary of V may be 
arbitrary and no useful coordinate system may be found, so other methods may need to 
be used. There are many different attempts. Amore pQ has applied a collocation method 
using so-called little sine functions for problems defined in two-dimensional domains of 
arbitrary shape. Chakraborty et al. [2] have presented an analytical perturbative method. 
In the two mentioned works brief surveys of other methods may be found. Recently, 
Steinbach et al. [3J have formulated a boundary element domain decomposition method 
that enables to transform the original problem to a new one defined on the boundaries 
separating the subdomains. 

The goal of this work is to present two methods that are applicable to the eigenproblem 
CO) - © in case the domain is such that it can be naturally divided into two non-overlapping 
subdomains. The methods consist in the application of a variational principle allowing 
the use of trial functions that may experience jumps in values or derivatives when passing 
from one subdomain to other. The Dirichlet-to-Neumann (DtN) integral operator or the 
Neumann-to-Dirichlet (NtD) integral operator is used, both are defined on the interface 
separating the subdomains. Each of the methods allows to replace the initial problem 
(PQ)-([2]) with a new problem defined in one of the subdomains and on the interface. The 
methods are related to the DtN and the NtD embedding methods for the bound states of 
the Schrodinger equation (defined in R 3 ) and their relativistic counterparts [HE]. The DtN 
method for the Schrodinger equation is a close relative of the embedding method proposed 
by Inglesfield [6j. In the Inglesfield's method the Green function formalism is used while 
in the DtN (and NtD) method an operator approach analogous to that employed in the 
R-matrix theory [?1[8] is applied. 

The structure of the paper is as follows. In section II a systematic construction of a 
variational principle (allowing the use of discontinuous trial functions) for bound states 
of the Helmholtz equation is presented. In section III the DtN and NtD operators are 
defined. Sections IV and V are devoted to the formalism of the DtN and NtD methods for 
bound states of the Helmholtz equation. In section VI a numerical example is provided. 

2 Variational principle allowing the use of discontinuous 
trial functions 

Let r be a two- (or more) dimensional finite domain of such a shape that it may be in a 
natural way divided into two subdomains, Tj and Tjj, separated by a smooth curve (or 
surface) denoted by S, as shown in figure CO Thus, the boundary of Ti consists of dTj 
and S while the boundary of Tu is composed of dTjj and S. A position vector lying on 
S will be denoted by p and n(p) will be the unit vector normal to S at the point p (we 
assume that n(p) is always pointed outward from Tj). Denoting 

*i(r) = *(r) (rer« ;i = /,Ji), (3) 

we may rewrite the initial problem ([I])-© as 

A*/(r) + fe 2 ^ 7 (r) = 0, reT/, (4) 

*/(r) = 0, r G dTj, (5) 
A*//(r) + k 2 ^ n (r) = 0, r G T n , (6) 
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Figure 1: Partitioning of the domain T into two subdomains Tj and Tjj, separated by the 
interface S; n(p) is the unit vector normal to the interface S at the point p. 



tf//(r) = 0, r G dV H . (7) 

The function \I/(r) and its gradient must be continuous in the whole domain T, so it is 
obvious that the functions ^/(r) and ^//(r) obey the equations: 

*j(p) - ttjj(p) = 0, (8) 

V±*/(p) - V ± *//(p) = 0, (9) 

where 

Vj.ti(p) = n(p) • V*i(r) (10) 

r=p 

is the normal derivative of at p. 

We want to determine the values of {k 2 } and the corresponding functions {^(r)}. 
Basing on equations ©, ©, © and © and using a method proposed by Gerjuoy et 
al. [9] we define a functional that provides some estimate of one of the sought values { k 2 } : 

F\k, W/j; Aj, A//, A, x] = k 2 + <A 7 |[A + F]*/), + <A W | [A + P]*//)^ 

+(A|*j- *//) + (x|V ± *j - V_l*jj). (11) 

The scalar products are defined as follows: 

($1$').= I dr$*(r)$'(r), (12) 

($|$') = J dp (p) (13) 

(dr is an infinitesimal volume element around the point r, dp is an infinitesimal scalar 
element of the interface S around the point p, and * denotes the complex conjugation). 
The value k, the function ^fj (vanishing on dTj) and the function ^fjj (vanishing on 
dTn) are some trial estimates of the exact quantities k, ^fi and Vl/77. The functions Aj 
(defined in Tj), A// (defined in Tjj), A and \ (both defined on S) play role of the Lagrange 
functions including equations (p j) , (j8 j) and ([9]) in the functional. The first variation 
of the functional (|lip with respect to arbitrary variations of k, VP/, about fc, ^j, Xl/j/ 
(supposing that the variations 5^>i and S^Sfu vanish on dTj and dTu, respectively) and 
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Aj, A//, A, x about some arbitrarily chosen A/, A//, A, \ may be written as 

£F[fc,¥/,tt//;A/,A Wj A,x] = 2fc<fA; [l + (A / |* / ) / + (A^)^}^] 
+ <[A + k 2 ]Aj\5^ I ) I + ([A + k 2 ]A II \5^ II ) n 
+ (A - V ± A/|5*/) - (A - V_lA//|5^//) 
+ (x + A/|V ± <5*/) - (x + Ajj|V±**//) 

+ (A J |V±<y*/) J + (A // |V ± ^//) JJ , (14) 

where 

($]$'). = A dr**(r)*'(r). (15) 

In the above scalar product, dr is an infinitesimal scalar element of dTi around the point 
r (cf. the definitions Q12|l and (|13|)). We seek such functions A/, Ajj, A and \ f° r which 
the functional is stationary, i.e. its first variation is equal to zero. So the functions A/, 
Ajj, A and x fulfil the equations: 

1 + (Aj|*j> i + <A 7/ |* /j ) // = 0, (16) 

[A + A; 2 ]Aj(r) = 0, r e r f , (17) 

Ai(r) = 0, re dT u (18) 

X(p)-V ± Ai(p) = 0, (19) 

X(p) + Ai(p) = 0, (20) 
where i = I, II From equations (119p and fj20f) we obtain 

Aj(p)-Ajj(p)=0, (21) 

V ± Aj(p) - V±Ajj(p) = 0. (22) 

Comparying equations (fI7|) . (flSjl. ([2T]) and ([22]) with equations flU)-© we find that A/(r) 
and A//(r) obey the same differential equations and the same boundary conditions as 
^j(r) and fyjj(r). This means that the functions Aj and Ajj are proportional to *$>i and 

Ai(r)=r/*i(r) (i = IJI). (23) 
The value of 77 may be found using the formulas (|23p in equation (I16p : 

7? = "<^ / |^ / > / + <^ // |* // > / / (24) 
According to equations (fTUj) . (|2U|) and (|23|) . we may write 

A(p) = r)[aV x Vl(p) + (1 - o)V±*//(p)], (25) 

X (p) = + (1 - 6)*jj(p)], (26) 

where a and 6 are arbitrary complex constants. 

Now, let us assume, that the trial Lagrange functions Aj, Ajj, A, x, appearing in the 
functional (jlip , are related to the estimates \P/ and VP// in the same way that the functions 
Aj, Ajj, A, x are related to the exact functions Vl/j and Vl///. Using the formulas obtained 
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from equations (J23D — d2SJ) by replacing the functions Vl/j, ^u, Aj, Ajj, A and x with the 
trial functions *$>i, A/, Ajj, A and x, transforms the functional (jlip to 







<W/|W J ) J H 


- (* w */y) 77 



(aV ± */ + [1 - a]V±*//|*/ - *//) 



yy 



(fefrj + [1 - 6]^jj|V ± ^j - V±^yj) 
+ <*/|^/> / + <*//|W // > // 

Our functional is supposed to estimate of a real quantity, so it should possess the property 

7*$ I ,^n]=Ff*i,^ii]- (28) 
After some rearrangements we find that equation (|28|) is obeyed if 

b = 1 - a* (29) 

and the final form of the functional is 

AW 7 ) j + (*/j|A* w ) w (aV±*j + [1 - a]V±*jj|*/ - *//) 



yy 



([l-a*]ttj + a*frjj|V_Lfrj-V ± frjj) 
(#/|*/) 7 + {^ji\^ii) h 

It is easy to verify that the exact functions are the stationary points of the functional (|30p : 

6F[* I ,* n ]=0, (31) 

and the corresponding stationary values are equal to k 2 : 

= k 2 . (32) 



The initial problem ([T])-([2]) is then equivalent to the variational principle (j3Q() — fj32|) . The 
important thing is that the functional (|30p allows to use trial functions \Pj and *$>u that, 
together with their gradients, are continuous in their domains, but do not have to match 
at 5, so at least one of the following cases may occur: 

+ *jj(p), Vj_*j(p) + V±*jj(p). (33) 

It is worth noting that such functionals are rather rarely applied. More details about the 
construction of similar functionals may be found in the paper of Szmytkowski et al. [10J, 
where variational principles for bound states of the Schrodinger and the Dirac equations 
have been presented. 

3 The DtN and NtD operators 

Let us assume that the subdomain Tu is such that we are able to find analytically the 
functions i/j(k, r) obeying the Helmholtz equation at some fixed real value of the parameter 
k (which need not be in the spectrum of the eigenproblem CE])-©): 

A^(«,r) + «V(«,r) = 0, reT n (34) 
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and the boundary condition 



ip(K,r)=0, r£dT n . (35) 

Now, let us define two mutually reciprocal integral operators B(n) and 7Z(k) such that for 
every i/}(k, r) at the interface S it holds that 

V ± iP(k,p) = B(k)iP(k, P ), (36) 

K(K)V±l/>(K,p)=l/>(K tP ) (37) 

(note that the operators Vj_ and B(k) are not identical, equation ([36]) is valid only for 
the functions if)(n,r)). The operator B(k) transforms the Dirichlet datum i/j(n,p) into 
the Neumann datum V ±ip(n, p) so it is called the Dirichlet-to-Neumann (DtN) operator. 
In analogy, the operator TZ(k) is called the Neumann-to-Dirichlet (NtD) operator. Using 
integral kernels B(k, p, p') and 7Z(k, p, p') of the operators, we may rewrite equations ([36]) 
and (j37|) in the following forms: 

V ± ^(k, P) = J s V B{k, p, pf)i/>(K, p'), (38) 

1&(k, P) = J °V K{K t p, f/)V±tl>(K, p'). (39) 
Now, let us analyze the eigensystem 

A^ n ( K ,r) + K 2 i(j n (K,r) =0, rer„, (40) 

^ n (K,r) = 0, redV u , (41) 

V_LVn(K, P) = b n (k)lj) n (K, p). (42) 

In the above system 6 n (re) is an eigenvalue, ip n (K, r) is an eigenfunction and k is some fixed 
real parameter. The eigensystem (j40]l -(j42l) is non-standard, because the eigenvalue 
appears not in the differential equation but in the boundary condition. Eigenproblems of 
such type are known as the Steklov eigenproblems 

Using the Green's theorem (and the condition (|4ip ) for two arbitrary eigenfunctions 
ip n (K,r) and ip n /(K,r) we obtain 

(ip n \Aip nl ) n - (AiJ; n \ip n/ ) n = (V ±ip n \ip n >) - (ipnlV ±ipn') ■ (43) 

In virtue of equation ()40j) the left-hand side of equation (|43p vanishes. Applying equation 
([HZ]) leads us to 

K(k)- 6„K«)](^«|^«') =0- (44) 

There are two conclusions we may draw from equation (1 44ft . First, if we take n' equal to 
n, we see that the eigenvalues are real 

b n (K) = b* n (K). (45) 

Second, if the eigenfunctions ip n (n,r) and ip n i(K,r) belong to different eigenvalues, they 
are orthogonal with respect to the surface scalar product (j 1 3 [) : 



(^n|Vw)=0 [b n (K) (46) 
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Now, let us assume that all the eigenfunctions of (|40|) ~(|42 p are orthonormal on S: 

(Vvi|Vv) = <W, (47) 

and that the surface functions {ip n (n, p)} form a complete set in the space of single- valued 
square-integrable functions defined on S and therefore obey the closure relation 

X}^(«,pX(«y) = 5 s (p-p'), (48) 

n 

where 5s(p — p') is the Dirac delta function on S. 

Combining the definition of the DtN operator (j36[) with equation (|42p we may write: 

B(K)tp n (K, p) = b n (K)i/) n (K, p). (49) 

We observe that eigenvalues of the DtN operator are the eigenvalues {6„(k)} of the eigen- 
system ([4"0jH(j4"2]) and the associated eigenfunctions are the surface parts {tp n (K, p)} of the 
eigenfunctions {tp n {n, r)}. According to equation (|38p . we may rewrite equation (|49p as 

J dp' B(k, p, p')ijj n (K, p') = b n (K)xl) n {K,p). (50) 

Multiplying the above formula by V'n(' s > p"), summing over n and using the closure relation 
(|48p leads us to the spectral expansion of the DtN operator kernel 

B(K,p,p') = ^2^ n (K, p)bn(K)lf>*(K, p'). (51) 
n 

As the DtN operator and the NtD operator are mutually reciprocal, the spectral expansion 
of the NtD kernel takes form 

K(k, p, p') = ^n(«, PK 1 ^*^, p')- (52) 

n 

It is obvious from the expansions (|51|) and (|52|) that the operators B(n) and TZ(k) are 
Hermitian. 

4 The DtN method 

If the trial functions employed in the functional (|30p are continuous on S, i.e. 

*/(p)=^(p), (53) 
the functional reduces to the following form: 

*W|»„¥„1 = JJA^hl^IlMlllu + W^-^I") . (54 ) 

Let us assume that the trial function ^jj is some function (re, r), obeying (|34"|) and 
CS]): _ 

* 77 ( r ) = ^( D )(K,r) (rer w ). (55) 
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Such choice of \F//(r) in virtue of equations (|36p and (|53p gives 

V ± *jj(p) = £(/c)*r(p). (56) 

Applying equations ([55]) . (plj) and ([56]) to equation leads us to such a form of the 
functional in which the only term containing iI)( d \k,y) is the integral (i/)( D ) \^ D ^) n '- 

Subtracting the complex conjugation of the Helmholtz equation for ^ d >{k,t) multiplied 
by dip( D '(K,r)/dK from the Helmholtz equation for ip( D >(K,r) differentiated with respect 
to k and multiplied by ip^ D '*(K,r) we obtain 

^ g )*(«,r)A ^ ( ° a )(/C>r) - ^ (D) (K,r) A ^ (D> = _ 2 ^p> ( } , (58) 

OK OK 



Integration of (|58p over I 1 // after employing the Green's theorem, the definition (|36p . the 
Hermiticity of 6 (ft) and the continuity constraint (|53p results in 



Substitution of equation (|59h transforms the functional (I57p into the form 



(59) 



<*i|*/>/ + &MraW*') 

which depends only on the parameter k and the trial function tyj defined in T/ and on 5. 
We arrive at a conclusion that the usage of the surface integral DtN operator allows us to 
reduce the initial problem defined in T to the subdomain Tj and the interface 5. 

We need to find such functions ^ d \k, r), (r £ T/), that make the functional (|60l) 
stationary. The associated stationary values are estimates of some of the values {& 2 } 
appearing in the problem (pQ)-(j2]). In practice, it may be impossible to find ^ d \k,t) 
analytically. In order to obtain some approximate solutions let us represent the trial 
function ^/(r) as a linear combination of some basis functions 0u(r), defined in Tj: 

M 

*/W = E«AW (reri). (61) 



Applying (JUT} to O yields 



_(fl) r _ t _, a f A^(K)a 

^ |M ' al = atA(^) W a ' (62) 



where a is an M-component column vector with elements {a^} and a T is its Hermitian 
adjoint, while A^'(k) and A^)(k) are M x M Hermitian matrices with elements 

a£?(«) = [-<^|a^> 7 + (^|v ± &, - + , (63) 
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Ag } («) = + (64) 

Let a^(/t) and a"( D N(n) be such particular vectors a and that make the functional 
(I62j) stationary with respect to variations in their components: 



6^ D) [K,a^( K ),a (D ^{ K )} = 0. (65) 
From equations (|65p and (|62p we arrive at the algebraic eigensystem 

A( d \k)z( d \k) = f( d \k)A( d \ K )z( d \k) (66) 
(and its Hermitian conjugate), where F^ d \k) is defined as 

f(°)(k) = ^ D) [«,a( D )t(«),a( D )(«)]. (67) 

The eigensystem (|66p has < M eigenvalues {-F 7 ^(«)} and corresponding eigenvec- 

tors {a~ *- d ^(k)}. These eigenvalues are second-order variational estimates of some among 
values {A; 2 } of the system (JT|) — (J2J) . In virtue of equation ([B"Tj) . the eigenvectors {a^ ( d \k)}, 
with the components {aj^ (k)}, give us functions: 

M 

¥<PW) = £a<#(K)*„(r) (r e T 7 ), (68) 
n=l 

which are first-order variational estimates of some of the eigenfunctions of the system 
CO) - © in the subdomain Tj. Now we may find the functions {ipj («,r)} which are the 
estimates of the eigenfunctions in Tjj . Let us expand them in the basis constitued by the 
eigenfunctions of the Steklov system (j4"0j) -(j4"2l): 

Vf>W) = J2 c m^M ( r e r ")- ( 69 ) 

n 

Letting the point r tend to the interface S, employing the orthonormality relation (1471) . 
and using the formula ([53]) , we obtain 

c^)(«) = (^„|*^)). (70) 

5 The NtD method 

In the previous section we started our reasoning with the matching condition (|53p for the 
trial functions used in the functional (|30p . Now, let us turn to another possibility and 
impose a weaker condition 

V_l*j(p) = VjJnip). (71) 
In this case, the functional (|30p simplifies to 

(^/|*/> / + (^//|*//> // (*/ 1 ttj^ + ^jj \Vii) n 



We assume that the trial function ^//(r) is some function obeying the set 

tfjj(r) = ^ N \n,r) (reTu). (73) 
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Investigation analogous to that leading to (|59p yields 

(^n n = -i(^\§^) (74) 



and the functional (17211 transforms to 



(*/|^> / -^(v ± ^|[a7e/^]v ± ^) 

Following the method of algebraization applied in case of the DtN method (see the formulas 
(I6ip - (l66p ) we arrive at the generalized matrix eigensystem 

AW( K )lW( K ) = F W(«)AW(«)a ^(«) (76) 

(and its Hermitian matrix conjug ate), where A^(k) and A( n \k) are M x M matrices 
with elements 

AgP(«) = [-<^|A^> 7 + (V ± ^\KV ± ^ -</> v - |[0&/0k]V_l&,)] (77) 

and 

Aj'W = (<^)/ - ^(V_ L AI |[0fc/0K]V_L0„). (78) 

Eigenvalues {F 7 ^(k)} of ([76]) are second-order variational estimates of some of the val- 
ues {k 2 } appearing in the set dH) — ([2]) , while components of the associated eigenvectors 
{a 7 ( n '(k)} yield the estimates of eigenfunctions of (P)-© in Tj: 

M 

*f )(«,r) = ^J»( K )^(r) (r G T z ). (79) 

The last step is to find estimates of eigenfunctions in Tjj. We expand i/j^(k, r) as follows 
^(«,r) =E c S ) («)^(«,r) (r G r I7 ). (80) 

n 

The orthonormality relation ()47|) . the properties of the NtD operator and the matching 
condition (JTTJ) lead to 

C W(K) = C(«)(^|V ± *f)). (81) 

It is worth noticing that in general the DtN method and the NtD method will give 
different estimates of the solutions of the initial system. 

More details about the DtN and NtD methods (for bound states of the Schrodinger 
equation and the Dirac equation in R 3 ) may be found in the works of Szmytkowski and 
Bielski 01151 . 



6 Numerical example 

To test the two methods, a few series of numerical calculations have been performed. A 
system in which V is a two-dimensional domain consisting of a semicircle of radius a joined 
to a rectangle of sides a and b, as depicted in figure has been examined. The first step 
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Figure 2: Geometry of the system used in the numerical illustration. 



is to decide which part of the whole domain T should be Tj and which one should be 
Tjj. The decision depends on the simplicity of the construction of the DtN and the NtD 
operators. The region in which it is easier to solve (|40[) — (j42j) should be taken as Tjj. In 
our example the subdomain Tj is the semicircle and the subdomain Tjj is the rectangle. 
It is not difficult to verify that the eigenvalues of (|40[) — (|42j) in this case are 



b n (n) 



-\Jk 2 — n 2 ir 2 /4a 2 cot(y / K 2 — n 2 7r 2 /4a 2 6) for k 2 > n 2 7r 2 /4a 2 
-y / n 2 7r 2 /4a 2 — n 2 coth(Y / re 2 7r 2 /4a 2 — k 2 b) for k 2 < n 2 it 2 /Aa 2 



(82) 



where n = 1, 2, The corresponding eigenfunctions are of the form 

ip n (K,x,y) = 

mr(x + a 



A n sin 



2a 




i/k 2 - n 2 7r 2 /4a 2 (y + b) 
\J n 2 ir 2 /Aa 2 — k 2 (y + b) 



for k 2 > n 2 TT 2 /Aa 



for k 2 < n 2 7r 2 /4a 



(,83) 



where A n according to the relation (|47|) are 

\J k 2 — n 2 7T 2 /Aa 2 b ^ 
^Jn 2 ^ 2 1 4a 2 — n 2 b 



/ a sin 
/~a sinh 



for k > n 2 n /4a 
for k 2 < n 2 7r 2 /4a 



(84) 



Using equations ([82])-([8l]) in and ([52}, we obtain the kernels of the DtN and NtD 
operators. Let us observe that in the examined case we may distinguish the even (sym- 
metric with respect to the y-axis) and the odd (antisymmetric) modes. We may search for 
them separately (which means working with smaller matrices), applying apropriate basis 
functions {^(r)} in (|6ip. For the even modes we may use 



0i(r) = r - a, r G T/, 

M (r) = 4> nm (r) = r sin[na(r - a)] cos(m/3</?), r G T/ 
(/i = 2,3,..., n,m = l,2,...) 

and for the odd modes we may apply 

M (r) = nm (r) = r sin[na(r - a)] sin(m/3(p), r £ Tj 
(/x=l,2,..., n,m=l,2,...) 



(85) 
(86) 



(87) 
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(each /i represents an unique combination of two integers n and m). In the above formulas 
ip is the angle between the y-axis and the position vector r = [x, y] (we assume that tp is 
positive for x < and negative for x > 0), r is the length of r, while a and j3 are some 
arbitrary real parameters. Note that all the functions vanish on dTj. The variational 
bases are formed from the functions with 1 < n < n max and 1 < m < m max . Do not 
forget about the extra function (|85p used for the symmetric states. The function is added 
to the basis because all the functions (|86p are equal to zero at r = 0, while in general the 
eigenfunctions of the even modes may be nonzero at r = 0. 

To obtain the estimates of some k 2 and Vl/(r), we must establish an initial value of k 
(which is some estimate of k). Then we apply some chosen variational basis, calculate 
the matrix elements of A^ d \k) and A^ d \k) or M n ^(k) and A^ n \k) and solve the matrix 
system ([66]) or (fT6j) . The resulting eigenvalues {F 7 ( D or n \k)} are used to set new values 
of k (separately for each 7). We focus ourselves on an arbitrary chosen state, let it be the 
state with 7 = 7', so we take 




We find new matrices, solve the new matrix system and obtain new estimates of eigenvalues 
{Fj^ DotN \k)} with Fy( DorN \K) among them. We then apply (JHSJ again and the 
iterative procedure repeats until convergence of Fy ( D or n \k) is achieved. 



Iteration 


k (D} 

even,l 


k (N) 

even,l 


^odd,l 


k (N) 


1 


2.0633 


2.0487 


3.4586 


3.4200 


2 


2.0611 


2.0604 


3.4508 


3.4447 


3 


2.0611 


2.0611 


3.4507 


3.4505 


4 




2.0611 


3.4507 


3.4507 


5 








3.4507 



Table 1: Convergence rate of the DtN and the NtD variational estimates of k of the 
lowest even mode and the lowest odd mode of the system used in the numerical example. 
The results obtained by employing the basis functions (|85p ~ (|87p with 1 < n < 15 and 
1 < m < 15. The inputs for the iteration procedure have been k = 2.0116 for the first 
even mode and k = 3.3836 for the first odd mode. SI units are used. 

The numerical calculations have been made for a = 1, b = 1.5 and a = /3 = 1 (SI units 
are used) . Table CD presents estimates of k of the first even mode and the first odd mode 
obtained by using the iterative procedure described above. The variational bases have 
contained the functions ([85]) and ([86]) (for the even mode) or ([87]) (for the odd mode) with 
1 < n < 15 and 1 < m < 15. The initial values of k have been k = 2.0116 for the first 
even mode and k = 3.3836 for the first odd mode, they are the exact values of k (with 
accuracy of 4 decimal places) of the first even mode and the first odd mode of the system 
([T])-([2]) with r being a rectangle of sides 2 and 2.5 (in the next series of calculations the 
fact that in case of the rectangle the second even mode is characterized by k = 2.9638 and 
the second odd mode by k = 4.0232 has been used). We see that the estimates converge 
after a few iterations, and that the results obtained by the DtN method converge faster. 
To verify, how the converged (after the iterative procedure) estimates of k depend on the 
basis dimensions, let us analyze the results collected in table [2j The two lowest even 
modes and the two lowest odd modes of our problem have been examined. The inputs 
for the iterative procedures have been k = 2.0116 and k = 2.9638 for the even modes and 
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W"max i^max 


even,l 
even,l 


i (N) 
^even,2 


"'odd,! 


^odd,2 
Vfd,2 


3,3 


2.0630 
2.0628 


3.0745 
3.0809 


3.4527 
3.4527 


4.2234 
4.2234 


5,5 


2.0611 
2.0611 


3.0734 
3.0734 


3.4511 
3.4511 


4.2200 
4.2200 


15,15 


2.0611 
2.0611 


3.0731 
3.0731 


3.4507 
3.4507 


4.2190 
4.2190 


25,25 


2.0611 
2.0611 


3.0730 
3.0730 


3.4506 
3.4506 


4.2189 
4.2189 


30,30 


2.0611 
2.0611 


3.0730 
3.0730 


3.4506 
3.4506 


4.2189 
4.2189 



Table 2: Converged DtN and NtD variational estimates of k of the two lowest even modes 
and the two lowest odd modes. The basis functions ([85|) - ([87|) with 1 < n < n max and 
1 < m < m max have been used. The inputs for the iterative procedures have been 
K = 2.0116 and k = 2.9638 for the even modes and k = 3.3836 and k = 4.0232 for the odd 
modes. SI units are used. 



k = 3.3836 and k = 4.0232 for the odd modes. It is seen that even quite small bases lead 
to estimates of good accuracy. The aim of the next series of numerical calculations was to 
estimate eigenfunctions $ of the four states. Figure [3] shows density plots of l^l 2 of the 
four modes obtained by using the DtN method in a way presented at the end of section 
IV. In the domain Tj the basis functions (I85p ~ (l87p with 1 < n < 25 and 1 < m < 25 have 
been used. The values of k have been equal to the converged estimates of k (see table [2]). 
Density plots obtained by the NtD method are exactly the same. 
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Figure 3: Density plots of |*| 2 : (a) and (b) - the two lowest even modes; (c) and (d) - 
the two lowest odd modes. The results obtained by the DtN method (see section IV). The 
basis functions (j85p -( j87]) with 1 < n < 25 and 1 < m < 25 have been used. The converged 
estimates of k (see the last row of table [2|) have been taken as the values of k. 
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